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We propose to parametrize the Stillinger-Weber potential for covalent materials starting from the 
valence force field model. All geometrical parameters in the Stillinger-Weber potential are deter¬ 
mined analytically according to the equilibrium condition for each individual potential term, while 
the energy parameters are derived from the valence force field model. This parametrization approach 
transfers the accuracy of the valence force field model to the Stillinger-Weber potential. Further¬ 
more, the resulting Stilliinger-Weber potential supports for stable molecular dynamics simulations, 
as each potential term is at energy minimum state separately at the equilibrium configuration. 

We employ this procedure to parametrize Stillinger-Weber potentials for the single-layer M 0 S 2 and 
black phosphorous. The obtained Stillinger-Weber potentials predict accurate phonon spectrum and 
mechanical behaviors. We also provide input scripts of these Stillinger-Weber potentials used by 
publicly available simulation packages including GULP and LAMMPS. 
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I. INTRODUCTION 

The atomic interaction is a fundamental ingredient for 
numerical investigation of nearly all physical or mechani¬ 
cal processes. For instance, in molecular dynamics (MD) 
simulations, the atomic interaction provides the retract¬ 
ing force for each atom in the Newton’s equation. There 
have been huge number of available potential models for 
the atomic interaction within different materials. For the 
covalent material, some representative potential models 
are shown in Fig. [1] in the order of their simulation cost; 
i.e., valence force field (VFF) model, Stillinger-Weber 
(SW) potential, Tersoff potential, Brenner potential, and 
ab initio approaches. These potentials (or approaches) 
are able to describe the bond stretching and angle bend¬ 
ing motions, which are two dominant motion styles in 
covalent materials. The bond twisting motion can also 
be treated by these potentials, although the twisting en¬ 
ergy is usually very small. 

The VFF model is a linear model, and is suitable for 
analytic derivation of many elastic quantities, so this 
model requires only limited computation cost. As an 
advantage of the VFF model, its parameters can be de¬ 
termined of high accuracy by fitting directly to some ob¬ 
servable elastic quantities. As a result, the VFF model 
was very popular for covalent materials, especially be¬ 
fore 1980s, when the CPU speed was very low. Conse¬ 
quently, the VFF model for most covalent materials have 
been well developed. For instance, the VFF model for 
M 0 S 2 has been proposed in 1975?i while the VFF model 
for black phosphorus (BP) was proposed in 1982;^ and 
the VFF model for graphene was developed in 1990 by 
Aizawa et alj^ These VFF models are useful for the study 
of many elastic properties in these quasi-two-dimensional 
nano-materials in recent years, especially during the gold 



FIG. 1: A schematic diagram comparing the simulation cost 
of different atomic interactions; i.e., VFF model, SW poten¬ 
tial, Tersoff potential, Brenner potential, and ab initio ap¬ 
proach. 


rush of graphene in the past decade. 

While the VFF model is beneficial for the fastest nu¬ 
merical simulation, its strong limitation is the absence 
of nonlinear effect. Due to this limitation, the VFF 
model is not applicable to nonlinear phenomena, for 
which other potential models with nonlinear components 
are required. The ab initio approach is accurate and ap¬ 
plicable to nonlinear phenomena, but it requires the most 
expensive simulation cost, due to the solution of the full 
quantum electronic problem. However, this approach de¬ 
sires the most expensive simulation resources. As a re¬ 
sult, the ab initio approach usually cannot simulate more 
than around a few thousand atoms, which poses serious 
limitations for comparisons to experimental studies. 

We are now aware that the VFF model is the cheapest 
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in computation cost, but it only works for elastic proper¬ 
ties. On the other hand, the ab initio approach can simu¬ 
late nearly all physical processes with high accuracy, but 
it requires the most expensive computation cost. Hence, 
the bridging between these two extreme cases is of prac¬ 
tical significance, since lots of studies prefer efficient sim¬ 
ulation with reasonable accuracy for the nonlinear treat¬ 
ment. There have been several potential forms to fill 
this bridging domain; including SW potentialr'^— Tersoff 
potential)2r2^ and Brenner potential All of these po¬ 
tential forms comprise reasonable accurate nonlinear ef¬ 
fects, and are particularly suitable for MD simulations. 

Among these potentials, the SW potential is one of the 
simplest potential forms with nonlinear effects includedi^ 
An advanced feature for the SW potential is that it in¬ 
cludes the nonlinear effect, and keeps the numerical sim¬ 
ulation at a very fast level. As a result, the SW potential 
has been widely used in the numerical simulation com¬ 
munity. The SW potential was originally proposed by 
Stillinger and Weber to describe the interaction in solid 
and liquid forms of silicon, and it has been used in other 
covalent materials like single-layer M 0 S 2 (SLMoS 2 )^ and 
single-layer BP (SLBP)i^ 

For chemically different materials, the SW potential 
form keeps unchanged, but all parameters need to be de¬ 
termined properly. In all present works, the parametriza- 
tion of SW potential (and also Brenner and Tersoff poten¬ 
tials) are done by fitting to some experimentally known 
quantities like the Young’s modulus, phonon spectrum, 
cohesion energy, and etc. Actually, from the above dis¬ 
cussion, we have learnt that most covalent materials al¬ 
ready have an accurate VFF model, which can describe 
linear properties accurately. Such attractive essence 
should be helpful for the parametrization of atomic po¬ 
tentials like SW potential, Tersoff potential, and Brenner 
potential. However, to-date, the accuracy of the VFF 
model was not transferred to other atomic potentials 
during their parametrization process. The present work 
takes the SW potential as an example to demonstrate 
the relationship between the VFF model and the SW 
potential. In doing so, we illustrate that the SW poten¬ 
tial parameters can be analytically parametrized based 
on the VFF model. 

In this paper, we propose a parametrization procedure 
for the development of SW potentials based on the VFF 
model. All SW geometrical parameters are determined 
according to the equilibrium condition for each SW term, 
while the SW energy parameters are derived from the 
VFF model analytically. This parametrization procedure 
is employed to develop the SW potentials for SLM 0 S 2 
and SLBP, which provide accurate phonon spectrum and 
mechanical behaviors. 

The present paper is organized as follows. In Sec.B, we 
present details about the parametrization of SW poten¬ 
tial based on the VFF model. The parametrization pro¬ 
cedure is applied to develop the SW potential for SLM 0 S 2 
in Sec.IB. Sec.IV is devoted to the analytic parametriza¬ 
tion of the SW potential for the SLBP. The paper ends 




FIG. 2: Two typical interactions in covalent materials. Each 
interaction term can be described using the VFF model or 
the SW potential, (a) Two-body bond stretching interaction, 
(b) Three-body angle bending interaction. Atom moving di¬ 
rections are depicted by red arrows. 

with a brief summary in Sec.V. 

II. VFF MODEL AND SW POTENTIAL 

For most covalent bonding materials, the bond stretch¬ 
ing and the angle bending are two typical motion styles 
as shown in Fig. [21 The corresponding interactions can 
be described by the VFF model in the linear regime for 
small bond variation Ar and angle variation A0, 


Vr = 

i/4 {Arf , 

(1) 

Vg = 

i/4did2 (Ad)2, 

(2) 


where and Kg are two VFF parameters. The W term 
is the potential that captures a variation in the bond 
length Ar. The Vg is for the potential corresponding to 
the variation of the angle Ad, where the anlge 9 is formed 
by two bonds of length di and ^ 2 . 

Besides VFF model, the SW potential is another use¬ 
ful potential for these two typical interactions in Fig. [2] 
There are two-body and three-body interactions in the 
SW potential, 

V 2 = (B/r^ - l) , (3) 

V3 = (cosd-cosdo)^ 

(4) 

where V 2 corresponds to the bond stretching and V 3 asso¬ 
ciates with the angle bending. The cut-offs Tmax, ?’maxi 2 
and rmaxi 3 are geometrically determined by the mate¬ 
rial’s structure. There are five unknown geometrical pa¬ 
rameters, i.e., p and B in the two-body V 2 term and pi, 
P 2 ^ and do in the three-body V 3 term, and two energy 
parameters A and K. 

Let’s assume that the material’s structure (bond length 
d and angle do) has been identified via experiments or 
other accurate theoretical methods. Using these knowl¬ 
edge, we can determine geometrical parameters in the 






3 


SW potential. First of all, it is reasonable to require that 
all bonds are at their equilibrium length and all angles 
are at their equilibrium angle value in the equilibrium 
configuration. That is, we have the equilibrium condi¬ 
tion, ^^\r=d = 0 and ^^\0=Bg = 0, for each bond and 
each angle individually. From ^^\r=d = 0 , we obtain the 
following constraint for parameters p and B in V2, 


-AB {d - rmaxf 
(Bd - rfS) 


( 5 ) 


where d is the equilibrium bond length from experiments. 
Hence, there is only one free geometrical parameter left 
in V2. In other words, Eq. JS]) ensures that the bond has 
an equilibrium length of d and the V2 interaction for this 
bond is at the energy minimum state at the equilibrium 
configuration. 

The three-body V3 term shown in Eq. o ensures 
= 0 explicitly, so we have no constraint on geomet¬ 
rical parameters for the three-body term. In fact, there 
is no free geometrical parameter in V3, because the an¬ 
gle 60 is from the experiment while pi and p2 have been 
determined by Eq. dSD- 

The energy parameters A and K in the SW potential 
can be derived from the VFF model, by equating the 
force constants from SW potential and the force con¬ 
stants in the VFF model. More specifically, we have 
^^\r=d = Kr and ^^\b=0o = K0did2 at the equilib¬ 
rium structure, leading to. 


A = 


K = 


Kr 

KBdid2 

2 Sin^ ^ ~^maxl2)-t"P2/ “^maxl3 )] 


where the coefficient a in Eq. m is, 


a = 


{d ‘f^max) 
P 

id T max) 


{B/d^ 
{B/d^ - 




-1) 

1 ) 



( 6 ) 

( 7 ) 


( 8 ) 


The bond length of the arms for the angle are c?i and c?2, 
which are from experiments or other theoretical calcula¬ 
tions. As a result, energy parameters in the SW potential 
are analytically related to the energy parameters in the 
VFF model. 

We summarize the key steps in the above analytic 
parametrization of the SW potential. In the SW poten¬ 
tial, bond stretching interaction is described by Eq. ©, 
and angle bending interaction is described by Eq. ®. 
The potential parameters are determined in three steps. 
First, interaction cut-offs (r^ax, ?’maxi2, and r^axis) are 
determined geometrically by the equilibrium configura¬ 
tion of the material. The bond length (d, di, and ^2) 



FIG. 3: (Color online) Atomic configuration of SLM 0 S 2 . 
There are two interaction types, i.e., the bond stretching term 
(red online) and the angle bending term (blue online). The 
x-axis is in the armchair direction, and the y-axis is in the 
zigzag direction. 


and the angle (ffg) are also from the experiment or other 
theoretical calculations. Second, geometrical parameters 
p in the two-body term and pi and p2 in the three-body 
term are determined by Eq. m, by assuming that each 
two-body SW term is at equilibrium separately. Third, 
energy parameters (A and K) are determined by Eqs. © 
and based on the VFF model. In this way, we have 
analytically determined nearly all SW potential parame¬ 
ters uniquely, except the parameter B for two-body SW 
potential in Eq. ©. The above derivation shows that 
there is no constraint imposed on the parameter B in 
the linear regime. The only condition for B to satisfy is 
that B < d'^, so that p > 0 . We will explain in the next 
two sections that the parameter B is related to the non¬ 
linear mechanical process, and should be fixed according 
to a nonlinear quantity. 


Before further processing, we note some advantages 
for the SW potential derived in this approach. First, 
such SW potential has fully inherited the accuracy of the 
VFF model, so it provides accurate description for linear 
properties which can be accurately described by the VFF 
model. Second, the equilibrium structure has been pre- 
built-in during the derivation as shown by Eq. , so 
this SW potential gives accurate relaxed configuration 
intrinsically. Third, each two-body and three-body term 
in the SW potential is fully relaxed separately at the 
equilibrium configuration; i.e., all bonds and angles are 
relaxed individually at the relaxed configuration. Hence, 
the SW potential will be extremely stable during MD 
simulations. Fourth, the SW potential includes nonlinear 
effects through the nonlinear forms of both two-body and 
three-body terms as shown in Eqs. m and 0, so the 
SW potential is able to provide nonlinear properties, eg. 
via performing MD simulations. 
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TABLE I: The VFF model parameters for SLM 0 S 2 from 
RefQ. 


Kr (4L) 


Kr, if) 

8.640 

0.937 

0.862 


TABLE II: Two-body (bond stretching) SW potential param¬ 
eters for SLM 0 S 2 used by GULP. The expression is V 2 = 
(B/r^ - l). 


A (eV) p B (A'*) 

^max (-^) 

Mo-S 

6.918 

1.252 17.771 0.0 

3.16 


III. SW POTENTIAL FOR MOS 2 

As an example, we apply the above parametrization 
procedure to develop the SW potential for SLM0S2 in this 
section. We use the equilibrium structure for SLM0S2 
from the first-principles calculations as shown in Fig. |31 
The bond length between neighboring Mo and S atoms is 
d = 2.382 A, and the angles are 9 = LSMoS = 80 . 581 ° 
and V' = /-MoSMo = 80 . 581 °. 

The VFF model for SLM0S2 is from Ref[l|, which is 
able to describe the phonon spectrum and the sound ve¬ 
locity accurately. We have listed the first three lead¬ 
ing force constants for SLM0S2 in Tab. U neglecting 
other weak interaction terms. The bond stretching term 
is 14 = -^ with Ad as the length variation of 

Mo-S bond (eg. Moi-Si). The angle bending term is 
Vg = ^d? (A 0 )^ for the angle Mo-S-S with Mo as the 

apex (eg. LSi^MoiSg,), and Vy, = for angle 

S-Mo-Mo with S as the apex (eg. /.MoiSqMos). 

Using Eqs. (0, ®, and 0, we obtain the SW poten¬ 
tial parameters for SLM0S2 used by GULP— as listed in 



FIG. 4: (Color online) Phonon spectrum for SLM 0 S 2 along 
the FM direction in the Brillouin zone. The results from the 
SW potential (lines) are compared with the experiment data 
(pentagons) from RefQI. The parameter B has no effect on 
the phonon spectrum. 



FIG. 5: (Color online) The effect of parameter B on the stress- 
strain relation for SLM 0 S 2 of dimension 27.0 x 28.1 A along 
the armchair direction at 1.0 K. The stress-strain curve is fit¬ 
ted to function a — Ee -|- , with E as the Young’s mod¬ 

ulus and D as the TOEC. The left top inset shows that the 
parameter B has no effect on the elastic property, Young’s 
modulus; while the right bottom inset shows that the pa¬ 
rameter B dominates the nonlinear quantity, TOEC, which 
is fitted by function D = —2953.8B^. The blue circle in the 
right bottom inset represents D = —899.8 GPa from the first- 
principles calculation,— which fixes parameter B = 0.552d^ 
for the SW potential. 



FIG. 6: (Color online) Stress-strain for SLM 0 S 2 of dimension 
27.0 X 28.1 A along the armchair and zigzag directions. The 
Young’s modulus is the same in the armchair and zigzag di¬ 
rections. The nonlinear mechanical properties are anisotropic 
in the armchair and zigzag directions. 

Tabs. Inland IIIII We have found in Sec.II that the param¬ 
eter B can not be determined by the linear VFF model, 
because B corresponds to the nonlinear mechanical be¬ 
havior. In other words, parameter B has no effect on 
linear properties. For instance, we compute the phonon 
spectrum for the SLM0S2 using two different sets of SW 
potential with B = O.ld^ and B = 0 . 552 d^. Although 
these two SW potential sets look completely different, 
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TABLE III: Three-body (angle bending) SW potential parameters for SLM0S2 used by GULP. The expression is V 3 = 
j^g[pi/(ri 2 -rmaxi 2 )+P 2 /(ri 3 -rmaxi 3 )l (cQS^ — cos^o)^- Mo-S-S iudicates tile bending energy for the angle with Mo as the apex. 



K (eV) 

6*0 (degree) 

Pi (A) 

P 2 (A) 

rminl 2 (A) 

^maxl2 (-^) 

rminlS (A) 

^maxl3 (-^) 

rmin23 (A) 

^max23 

Mo-S-S 

67.883 

81.788 

1.252 

1.252 

0.0 

3.16 

0.0 

3.16 

0.0 

3.78 

S-Mo-Mo 

62.449 

81.788 

1.252 

1.252 

0.0 

3.16 

0.0 

3.16 

0.0 

4.27 


TABLE IV: SW potential parameters for 
pression is V2 = eA e 

(^^S 0,.,^ - COS 0o)^ 

LAMMPS. 


SLM0S2 used by LAMMPS.— The two-body potential ex- 
J . The three-body potential expression is V3 = 
The quantity tol in the last column is a controlling parameter in 



e(eV) 

G (A) 

a 

A 

7 

cos 60 

A 

Bl 

P 

q 

tol 

Mo-S-S 

1.000 

1.252 

2.523 

67.883 

1.000 

0.143 

6.918 

7.223 

4 

0 

0.0 

S-Mo-Mo 

1.000 

1.252 

2.523 

62.449 

1.000 

0.143 

6.918 

7.223 

4 

0 

0.0 


Fig. 0] shows that the phonon spectrum corresponding to 
different parameter B are exactly the same. 

To fix parameter B, a nonlinear quantity is needed. 
Fig. [5] clearly demonstrates that the parameter B has 
strong effect on the nonlinear mechanical behavior of the 
stress-strain relation during the tension of a SLM 0 S 2 of 
dimension 27.0 x 28.1 A at 1.0 K. The stress (cr) is fitted 
as a function of strain (e), cr = Ee + ^De^, with E as 
the Young’s modulus and D as the third-order elastic 
constant (TOEC). The left top inset in Fig.[5]shows that 
the parameter B has no effect on another elastic property, 
the Young’s modulus. Fig. [5] right bottom inset shows 
the relationship between D and parameter B. Using the 
first-principles result,— D = —899.8 GPa, we can fix the 
parameter B = 0.552d^. 

The SW potential parameters for SLM 0 S 2 used by 
LAMMPSi^ are listed in Tab. IIVI The potential script for 
LAMMPS can be found in the supplemental material— 
We use LAMMPS to perform MD simulations for the 
mechanical behavior of the SLM 0 S 2 under uniaxial ten¬ 
sion at 1.0 K and 300.0 K. Fig. [5] shows the stress- 
strain curve during the tension of a SLM 0 S 2 of dimen¬ 
sion 27.0 X 28.1 A. Periodic boundary conditions are ap¬ 
plied in both armchair and zigzag directions. The struc¬ 
ture is thermalized to the thermal steady state with the 
NPT (constant particle number, constant pressure, and 
constant temperature) ensemble for 100 ps by the Nose- 
Hoove r^^'^^ approach. After thermalization, the M 0 S 2 
is stretched in one direction at a strain rate of 10® s“^, 
while the stress in the lateral direction is allowed to be 
relaxed to be zero. We have used the inter-layer space in 
bulk M 0 S 2 , 6.092 A, as the thickness of the SLM 0 S 2 in 
the computation of the strain energy density. 

In Fig. ini from the curve in the linear region, e G 
[0, 0.01], we get the Young’s modulus of SLM 0 S 2 around 
165.7 GPa and 167.0 GPa in the armchair and zigzag 
directions, respectively. The shear modulus and Pois¬ 
son’s ratio can also be obtained in this linear regime. 
It is obvious that the Young’s modulus is isotropic for 


TABLE V: The VFF model parameters for SLBP from RefQ. 


Kr (^) 

Ke (4L) 

(^) 

7.578 

0.818 

0.710 


SLM 0 S 2 due to the three-fold rotational symmetry in 
this quasi hexagonal lattice structure— Recent experi¬ 
ments have measured the effective Young’s modulus to 
he E = 120 ± 30 or A = 180 ± 60 Nm'^,^ 

These values correspond to an in-plane Young’s modu¬ 
lus of 198.6 ± 49.7 GPa or 297.9 ± 99.3 GPa, considering 
the thickness of 6.092 A. Our theoretical values are quite 
close to the first experiment. The TOEG in the zigzag 
direction is larger than that in the armchair direction, 
which agrees with the first-principles calculations— The 
SLM 0 S 2 yields at smaller strain at 300 K than 1.0 K for 
both armchair and zigzag directions. 

In 2013, the author has parametrized with collab¬ 
orators a SW potential set (SW 2 OI 3 -M 0 S 2 ) for the 
SLM 0 S 2 by fitting parameters to the experimental 
phonon spectrum<^ The present SW potential (SW2015- 
M 0 S 2 ) has fewer interaction components than the 
SW 2 OI 3 -M 0 S 2 potential. However, the phonon spec¬ 
trum from SW 2 OI 5 -M 0 S 2 potential can be as accurate 
as the SW 2 OI 3 -M 0 S 2 potential, because the present 
parametrization procedure transfers the accuracy of the 
VFF model to the SW potential. Furthermore, each in¬ 
teraction component in the present SW2015-MoS2 po¬ 
tential is at equilibrium invidually, which is more strict 
than the SW2013-MoS2 potential, in which the equilib¬ 
rium condition is satished overall among all interaction 
components. As a result, the SW2015-MoS2 potential is 
more stable for MD simulations. 
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FIG. 7: (Color online) Configuration of SLBP. Atoms are 
divided into the top group (atoms 1, 2, and 3) and the bottom 
group (atoms 4, 5, and 6). There are two interaction terms, 
the bond stretching term (red online) and the angle bending 
term (blue online). The x-axis is along the armchair direction, 
and the y-axis is along the zigzag direction. 


TABLE VI: Two-body (bond stretching) SW potential pa¬ 
rameters for SLBP used by GULP. The expression is V 2 = 

^g[p/C-r^ax)] ^5/^4 _ 


A (eV) p (A) B (A'*) r^in (A) 


(A) 


P-P 3.626 0.809 14.287 


0.0 


2.79 



FIG. 8: (Color online) The effect of parameter B on the stress- 
strain relation for SLBP along the armchair direction at 1.0 K. 
The stress-strain curve is fitted to function a — Ee -f \De^, 
with E as the Young’s modulus and D as the TOEC. Left 
top inset shows that parameter B has no effect on the elastic 
quantity. Young’s modulus. However, the right bottom inset 
shows that the parameter B has strong effect on the nonlinear 
property, TOEC, which is fitted to function D = —13.8 — 
227.15^. The blue circle in the right bottom inset represents 
D = —91.3 GPa from the first-principles calculation,— which 
helps to fix parameter B = 0.584d^ for the SW potential. 


IV. SW POTENTIAL FOR SLBP 


As another example, we apply the parametrization pro¬ 
cedure to develop the SW potential for SLBP in this sec¬ 
tion. The structure for SLBP shown in Fig. [3 has been 
identified by experiment P atoms are divided into the 
top group (including atoms 1, 2, and 3) and the bot¬ 
tom group (including atoms 4, 5, and 6 ). There are two 
bond lengths, i.e., the intra-group bond (eg. bond 1 - 2 ) 
di = 2.224 A and the inter-group bond (eg. bond 1-4) 
d 2 = 2.244 A. These two bond lengths are very close to 
each other, so it can be assumed that both bonds have 
the same length o^ d = 2.224 A. The intra-group angle 
(eg. Z213) is 9 = 96.359° and the inter-group angle (eg. 
Z314) is iIj = 102.09°. 

Tab. IVl lists the VFF model parameters for SLBP from 
RefU. The bond stretching potential between two neigh¬ 
boring P atoms isVr = - (Ad)^. We note that the intra¬ 

group bond and the inter-group bond essentially have the 
same stretching parameter)^ As a result, there is only one 
VFF model parameter for bond stretching potential. The 
angle bending potential is Vg = {A9y for the intra¬ 
group angle, and = ^d?' (A'f/')^ for the inter-group 
angle. These three terms make dominant contribution to 
the interaction for the SLBP, while other weak interac¬ 
tion terms have been omitted in the present work. As 
a compensate, these parameters in Tab. are different 
from the original value by an overall factor of 0.76. 

Using Eqs. ©, ©, and o, we obtain the SW po¬ 
tential parameters for SLBP used by GULP— as shown 
in Tabs. EH and Em The determination of B is illus- 



FIG. 9: (Color online) Phonon spectrum for SLBP along FM 
from the SW potential is compared to the data from the ah 
initio calculation.—. 


trated in Fig. [ 8 l The parameter B has no effect on the 
elastic property, the Young’s modulus, as shown by the 
left top inset in Fig. [SI However, the parameter B has 
strong effect on the nonlinear quantity, TOEC, which can 
be fitted to the function D = —13.8 — 227.15^. Using 
this relationship between the TOEC and parameter B, 
we obtain the parameter B = 0.584d^ corresponding to 
D = —91.3 GPa from the hrst-principles calculations;^ 
We note that D = —13.8 A 0 even for H = 0, as shown 
in the right bottom inset of Fig. jS] For B = 0, the only 
nonzero SW potential term is V 3 = A(cos0 — cos do)^, so 
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TABLE VII: Three-body (angle bending) SW potential parameters for SLBP used by GULP. The expression is V 3 = 
^g[pi/(''i 2 ->-miixi 2 )-l-P 2 /(>-i 3 -’'maxi 3 )] liues are for intra-group angles. The last two lines are 

for inter-group angles. 



K (eV) 

60 (degree) 

pi (A) 

P2 (A) 

rminl2 (A) 

^maxl2 (-'^) 

rminl3 (A) 

^maxl3 (-^) 

rmin23 (A) 

^max23 (-'^) 

Pt-Pt-Pt 

35.701 

96.359 

0.809 

0.809 

0.0 

2.79 

0.0 

2.79 

0.0 

3.89 

Pb-Pb-Pb 

35.701 

96.359 

0.809 

0.809 

0.0 

2.79 

0.0 

2.79 

0.0 

3.89 

Pt-Pt-Pb 

32.006 

102.094 

0.809 

0.809 

0.0 

2.79 

0.0 

2.79 

0.0 

3.89 

Pb-Pb-Pt 

32.006 

102.094 

0.809 

0.809 

0.0 

2.79 

0.0 

2.79 

0.0 

3.89 


TABLE VIII: SW potential parameters for SLBP used by LAMMPS. The two-body potential expres¬ 
sion is V 2 = eA e ]. The three-body potential expression is V 3 = 

eAe+T'°'(’"3fe ] (cos— cos0o)^- The quantity tol in the last column is a controlling parameter in 

LAMMPS. Pt indicates atoms from the top group, while Pb represents atoms in the bottom group. 



6(eV) 

CT (A) 

a 

A 

7 

cos Oo 

A 

Bl 

P 

<7 

tol 

Pt-Pt-Pt 

1.000 

0.809 

3.449 

35.701 

1.000 

-0.111 

3.626 

33.371 

4 

0 

0.0 

Pb-Pb-Pb 

1.000 

0.809 

3.449 

35.701 

1.000 

-0.111 

3.626 

33.371 

4 

0 

0.0 

Pt-Pt-Pb 

1.000 

0.809 

3.449 
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FIG. 10: (Golor online) Stress-strain for SLBP during ten¬ 
sion process. Highly anisotropic mechanical behaviors are 
observed in the armchair and zigzag directions. 

the nonzero residue, D = —13.8 GPa, originates from the 
nonlinear effect purely contributed by the angle bending 
interaction. This is different from SLM 0 S 2 results shown 
in the right bottom inset in Fig. [51 where D = 0 at B = 0. 
This difference can be attributed to the different space 
groups for SLBP {C 2 h) and SLM 0 S 2 (Ds^i). As a re¬ 
striction of the three-fold symmetry in the SLM 0 S 2 , the 
overall nonlinear effect from the angle bending vanishes. 

The phonon spectrum for the SLBP from the SW po¬ 
tential is shown in Fig. | 9 l The results from SW potential 
agrees quite well with the first-principles calculations 

SW potential parameters for SLBP used by 
LAMMPSi^ are listed in Tab. IVIIII The potential 
script for LAMMPS can be found in the supplemental 


materialj^ We use LAMMPS to perform MD simula¬ 
tions for the tensile behavior for the SLBP of dimension 
26.3 X 29.8 A at 1.0 K and 300.0 K. Fig. (TUI shows the 
stress-strain curves during the tensile deformation of 
the SLBP along the armchair direction and the zigzag 
direction. Periodic boundary conditions are applied 
in both armchair and zigzag directions. The structure 
is thermalized to the thermal steady state with the 
NPT (constant particle number, constant pressure, 
and constant temperature) ensemble for 100 ps by the 
Nose-Hooveri2ii^ approach. After thermalization, the 
SLBP is stretched in one direction at a strain rate of 
10 ® s“^, and the stress in the lateral direction is allowed 
to be fully relaxed. We have used the inter-layer space of 
5.24 A as the thickness of the SLBP in the computation 
of the strain energy density. 

In Fig. 1101 from the stress-strain curve in the strain 
range [0, 0.01], we obtain the Young’s modulus 33.5 GPa 
and 105.5 GPa in the armchair and zigzag directions, 
respectively. These values are close to the previously re¬ 
ported ab initio results, eg. 28.9 Nm“^ in the armchair 
direction and 101.6 Nm“^ in the zigzag direction from 
Refl^. The SLBP yields at smaller strain at 300 K than 
1.0 K for both armchair and zigzag directions. 

In a recent work, the author has parametrized with 
collaborators a SW potential set (SW2013-BP) for the 
SLBP by fitting parameters to the phonon spectrum 
from ab initio calculationsi^ The present SW potential 
(SW2015-BP) has fewer interaction components than the 
SW2013-BP potential. However, the phonon spectrum 
from SW2015-BP potential can be as accurate as the 
SW2013-BP potential, because the present parametriza- 
tion procedure transfers the accuracy of the VFF model 
to the SW potential. Furthermore, each interaction com- 

















































ponent in the present SW2015-BP potential is at equilib¬ 
rium invidually, which is more strict than the SW2013- 
BP potential, in which the equilibrium condition is sat¬ 
isfied overall among all interaction components. As a 
result, the SW2015-BP potential is more stable for MD 
simulations. 

As a final note, this work proposes a method to develop 
the SW potential based on the VFF model, and applies 
this parametrization approach to SLM 0 S 2 and SLBP. 
The parametrization procedure, represented in Sec.II, is 
actually applicable to the development of other atomic 
potentials for a wide range of covalent materials. It is 
quite obvious that the SW potential for other covalent 
materials can also be developed analogously. 

An important technical note. For the simula¬ 
tion of SLM0S2 by LAMMPS, one needs to re¬ 
compile the LAMMPS package with our modi¬ 
fied source file, pair_sw.cpp, in the supplemental 
materialj^ This helps to exclude angle bending 
for angles like /.S1M01S4 in Fig. [ 3 |, which is not 
considered in the present work. However, for the 
simulation of SLBP using LAMMPS, one must 
use the original LAMMPS package; i.e., use the 
original source file, pairsw.cpp. 


V. CONCLUSION 


In conclusion, we have proposed an approach to deter¬ 
mine the SW potential parameters based on the valence 
force field model. The SW potential developed following 
this approach inherits the accuracy of the VFF model 
in the description of linear physical properties. Further¬ 
more, the accurate equilibrium structure information is 
pre-built-in, and this potential is very suitable for stable 
MD simulations. Finally, the SW potential can be easily 
used in many available MD simulation packages such as 
GULP and LAMMPS. As two examples, we apply this 
parametrization technique to develop the SW potential 
for SLM 0 S 2 and SLBP, which are found to provide accu¬ 
rate phonon spectrum and mechanical properties. 
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